I 


m 



LOAN COPY: RETURN TO 
AFWL (WLIL-2) 
KIRTLAND AFB, N MEX 


CHEMICAL EQUILIBRIUM 
OF ABLATION MATERIALS 
INCLUDING CONDENSED SPECIES 


by C. W. Stroud and Kay L. Brinkley 

Langley Research Center 
Langley Station , Hampton, Va. 


NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 


WASHINGTON, D. C. 


AUGUST 1969 


tech library kafb, nm 



TECH LIBRARY KAFB, NM 



1. Report No. 2. Government Accession No. 

NASA TN D-5391 

3. Recipient's Catalog No. 

4. Title and Subtitle 

CHEMICAL EQUILIBRIUM OF ABLATION MATERIALS 
INCLUDING CONDENSED SPECIES 

5. Report Date 

August 1969 

6. Performing Organization Code 

7. Author(s) 

C. W. Stroud and Kay L. Brinkley 

8. Performing Organization Report No. 

L-5957 

9. Performing Organization Name and Address 

NASA Langley Research Center 
Hampton, Va. 23365 

10. Work Unit No. 

124-07-13-01-23 

1 1 . Contract or Grant No. 

13. Type of Report and Period Covered 

Technical Note 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D.C. 20546 

14. Sponsoring Agency Code 

15. Supplementary Notes 

16. Abstract 


A computer program for calculating the multiphase chemical equilibrium of an arbitrary system is 
described and the complete program is included as an appendix. This program has been used to calculate 
the chemical equilibrium of low-density phenolic nylon over a wide range of conditions. The effects of 
including carbon in the condensed phase on equilibrium composition and ablative performance are dis- 
cussed. Methods for incorporating the important effects of chemical equilibrium into ablation analyses 
are described. 


17. Key Words Suggested by Author(s) 
Ablation materials 
Chemical equilibrium 
Phenolic nylon 


18. Distribution Statement 


Unclassified - Unlimited 


19. Security Classif. (of this report) 


20. Security Classif. (of this page) 


21. No. of Pages 


22 . 


Price* 


Unclassified 


Unclassified 


41 


$3.00 


*For sale by the Clearinghouse for Federal Scientific and Technical Information 
Springfield, Virginia 22151 



CHEMICAL EQUILIBRIUM OF ABLATION MATERIALS 
INCLUDING CONDENSED SPECIES 

By C. W. Stroud and Kay L. Brinkley 
Langley Research Center 

SUMMARY 

A method for computing the chemical equilibrium of complex systems including con- 
densed species is described. The computer program developed to perform calculations 
based on this method is included as an appendix. This program is used to compute the 
chemical equilibrium composition of pyrolysis products of phenolic nylon over a tempera- 
ture range from 500° K to 3500° K. The inclusion of condensed species as possible equi- 
librium products yielded a composition that had a simple form below 2700° K that can be 
applied to many ablation materials. The implications that these equilibrium results have 
on ablation analysis are discussed and ways to include the significant parts of these results 
in approximate ablation analyses are outlined. 

INTRODUCTION 

The gases that result from the pyrolysis of a charring ablator have a beneficial 
effect on ablator performance by influencing surface recession, by blocking energy at the 
external surface, and by absorbing energy in the char layer. In theory, a consideration 
of the kinetics of the reactions of the pyrolysis gases as they pass through the hot char to 
the external surface would result in the precise prediction of the chemical state of the 
gases and the most realistic appraisal of the influence of the gases on ablator perfor- 
mance. However, the chemical kinetics of the pyrolysis gases are not generally con- 
sidered in detail because of the lack of information on reaction rates and the extensive 
computations that would be required. 

Alternatives to a detailed study of the reaction kinetics of the gases are to assume 
that the gas composition is frozen in the initial pyrolysis state or to assume that the gases 
react to the equilibrium state corresponding to each temperature encountered during pas- 
sage through the char layer. The frozen composition approximation is the simplest 
approximation, but leads to a conservative estimate of the ablator performance. On the 
other hand, the assumption of chemical equilibrium avoids the conservatism of frozen 
composition while providing a practical alternative to the detail of the kinetics case. How- 
ever, previous calculations of the chemical equilibrium of pyrolysis gases have generally 



included only gaseous species and have neglected condensed species such as carbon. (For 
example, see ref. 1.) 

The paper presents a computer program that was developed to include both gaseous 
and condensed species in the determination of equilibrium compositions by minimizing 
the Gibb's free energy. In addition, a technique for minimizing computer time and an 
approximate model based on the use of a few of the many possible species are described. 

Example results are given from calculations for the pyrolysis gases of phenolic 
nylon for the temperature range of 500° K to 3500° K and pressures of 0.1, 1, and 10 atmo- 
spheres. The inclusion of the condensed carbon species is shown to have a significant 
effect on the surface- recession and heat- accommodating ability of a charring ablator. 

SYMBOLS 

The units used for the physical quantities defined in this section are given in cgs 
units and in the International System of Units (SI) (ref. 2). 

a elemental composition matrix 

b elemental total in mixture, moles 

Cp effective specific heat, j/mole 

f free energy of individual species, J/mole 

F Gibb's free energy, J/mole 

H c enthalpy plus chemical energy, J/mole 

H enthalpy, J/mole 

i,j,k,£ integers 

Kp equilibrium constant 

m number of elements in mixture 

N number of gaseous species in mixture 

n mole fraction 
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p partial pressure, atmospheres 

Pj. total pressure, atmospheres 

Q number of condensed species in mixture 

R gas constant, 1.987 J/mole-°K 

T temperature, °K 

x number of moles 

y value of x at previous iteration 

a degree of dissociation 

/3 convergence criterion 

p density, kg/m3 

Superscripts: 

g denotes gaseous species 

s denotes condensed species 

EQUILIBRIUM COMPOSITIONS BY MINIMIZATION OF FREE ENERGY 

In a closed, chemically reacting system at constant pressure and temperature, the 
potential function which governs chemical changes is the Gibbs free energy. The condi- 
tion for chemical equilibrium is that the free energy must be a minimum. (See ref. 3.) 
The technique used here to minimize the free energy was developed by modifying the 
well-known method of reference 4 to include condensed species. This modification is 
similar to that given in reference 5. 

Formulation of Problem 

The Gibbs free energy of a mixture of N gaseous species and Q condensed spe- 
cies can be written as: 
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where ( is the molal standard free energy. 
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The free energy of a condensed species is 

f.s = x s(F\* 
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(3) 


The equilibrium condition can be obtained by minimizing the free energy of the mix- 
ture (eq. (1)) subject to the mass balance constraints. These constraints are given by the 
following equations: 

N Q 

^ a i£ x i g + ^ =b Z (Z= i’ 2 ’ • • • m ) ( 4 ) 

i=l 3=1 


where m is the number of elements in the mixture. 


Method of Solution 

At a given temperature and total pressure, the free energies of all possible mix- 
tures of a given set of species form a surface in the hyperspace the coordinates of which 
are the species. However, in practical cases some areas on this surface are not accessi- 
ble because the mole fraction that a particular species may take on is constrained by the 
total amount of each element available for the mixture and by the fact that negative mole 
fractions are not realizable. Thus the objective is to find that location on the accessible 
areas where the free energy is a minimum. The set of coordinates, or mole fractions, 
that define the minimum point constitute the desired equilibrium composition. 
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The method used here to find the point of minimum free energy is the method of 
steepest descent (ref. 6) applied to a quadratic representation of the free-energy surface. 

The solution is initiated by selecting an arbitrary set of mole fractions from which 
a point on the free-energy surface is computed. An n-dimensional parabola (where, n 
equals the number of species) is fitted to the surface at the point corresponding to the ini- 
tial selection of mole fractions and in the direction of steepest descent. The minimum of 
this parabola is obtained by using Lagrangian multipliers to deal with the mass balance 
constraints. The coordinates of the minimum, or the point on the parabola closest to the 
minimum but yielding a set of positive mole fractions, provides a new set of mole frac- 
tions from which a new point on the energy surface is calculated. A parabola is generated 
at the new point and the procedure is repeated until the change in the set of mole fractions 
is less than a prescribed value. 

If the constraints are not preventing the minimum of the parabola from being 
reached, this procedure will quickly converge to the neighborhood of the minimum free- 
energy point. In this paper, convergence was judged to be sufficient if the sum of the 
absolute values of the changes in mole fractions from one iteration to another was less 
than a very small number; for example, 10“ 5. 

If the restriction to a positive set of mole numbers is causing the mole fractions to 
change very slowly, the solution can be expedited by removing the species that creates the 
constraint and continuing the calculations on a new free energy surface. When the mini- 
mum free energy point is reached, the species that was eliminated may again be included 
to determine whether it exists in the equilibrium mixture. 

The question of uniqueness of solution has been considered in reference 7 where it 
is proven that all minima of the Gibbs free-energy function of multiphase mixtures are 
global minima. 

The inherent advantage of this method for determining equilibrium compositions is 
that no prior knowledge of the reactions that take place in a mixture is required. As a 
matter of fact, experience has shown that a poor initial selection of the mole fractions 
does not change the final solution nor add greatly to the solution time. The method also 
affords the flexibility of changing the initial set of species as desired. 

Computer Program 

The computer program for the previously described solution was written in Fortran 
6000 language and a listing of the program is included as appendix A. This listing is com- 
plete and contains all the subroutines necessary to run the program. Numerous comment 
cards are included in the listing to identify the various sections of the program. As com- 
piled, the program can accept 90 gaseous species and 10 condensed species at one time. 
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These species may be composed of up to 10 elements. However, the dimension statements 
could be changed so that a larger number of species and elements could be accommodated. 

The inputs required by the program are listed in appendix B. The program is, in 
general, self-contained. Once the thermodynamic data are obtained for all species, the 
only information required is related to the temperature and pressure involved and an ele- 
mental analysis of the mixture. In addition, the desired accuracy of iteration /3 must 
be specified. The iteration cycle is terminated when the sum of the absolute value of the 
changes in mole fraction is less than /3 from one iteration cycle to the next. A typical 
value used for /3 is 10“ 

The equilibrium mole fraction and weight fraction of each species are printed as 
output. The standard output also includes the enthalpy of the mixture as inserted initially 
and at equilibrium conditions. The sum of the chemical energy and enthalpy is also cal- 
culated along with the molecular weight of the mixture. The enthalpy calculated is the 
sum of the partial enthalpies of the gaseous components. Thus, 

H = £ n k H k (5) 

k=l 

Appendix C lists the species commonly found in ablation materials and sources of 
data on these materials. At present, there are 53 species for which reliable data have 
been found. Any of these species or new species can be added to a case merely by adding 
data cards to the input deck. The program will then include this species in the calcula- 
tions and identify it whenever it appears. 

The program consists of a framework of programed logic that has been built around 
the basic method of solution described in the previous section. This logical sequence 
makes possible the inclusion of all possible species (condensed and gaseous) in the initial 
mixture. The program automatically eliminates those species that exist in negligible 
quantities or attempt to take on negative mole fractions. The elimination process assures 
rapid convergence by removing the artificial constraint placed on the system initially by 
including those species in the mixture. Once the equilibrium condition is determined with 
the remaining species, the eliminated gaseous species are put back in the mixture. Thus, 
no species actually present in the mixture can be eliminated. 

The efficiency of the elimination process is typified by a calculation that required 
6000 iterations before elimination, but only 10 iterations after the elimination process 
was adopted. With the modifications, as included, the program is capable of reasonable 
computing speed. The chief advantage of the program is its flexibility and the ease with 
which new species can be added. The speed of solution is more a function of the number 
of elements and condensed species involved than the number of gaseous species. For 
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example, many ablation materials are limited to the four elements: carbon, hydrogen, 
oxygen, and nitrogen. 

Since the number of equations equals the number of elements plus the number of 
condensed species plus one, there are only six simultaneous equations to solve at each 
iteration, and a large number of possible species may be included without undue expendi- 
ture of computing time. 

Routine calculations of the type reported herein with 4 elements, 1 condensed specie, 
and approximately 40 gaseous species required about 1 second of computing time per tem- 
perature. Simpler systems, such as the check cases to be discussed later, require 300 to 
400 milliseconds of computer time to complete a calculation. Since there is no great 
computing time penalty for including a large number of species, everything for which data 
are available can be included initially. After some experience is gained, those species 
that are never significant can be removed from subsequent calculations simply by 
removing data cards from the input deck. 

Program Check Cases 

In order to illustrate the accuracy of the computer program, results are compared 
with simple mixtures for which exact equilibrium compositions can be determined or 
experimental data exist. First, the decomposition of methane into solid carbon and hydro- 
gen is considered. The equilibrium composition of these three species was calculated at 
1273° K with a total pressure of 7.1 atmospheres by using the computer program. The 
results indicate the following set of mole numbers: 0.1095 CH4, 1.7810 H2, and 
0.8905 C(s) where (s) denotes the solid state. 

The exact equilibrium composition may be computed by using the equilibrium con- 
stant Kp. Consider the reaction 

CH 4 - C(s) + 2H 2 


for which by definition 



( 6 ) 


For ideal gases, the mole fraction 
therefore, 


n k equals the partial pressure traction P^P,; 


K " H * 2Pt 
P n CH, 


(7) 
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The mixture of species that results from the dissociation of methane can be symbolized 
as follows: 


CH 4 - QiC(s) + 20-H2 + (1 - a-)CH 4 (8) 

From this relation, it can be seen that there are 1 + a moles of gas species and that 

2a 




n n 2 = r 


+ a 




n 


CH 


1 - a 

4 \ + a 


(9) 


If equations (9) are substituted into equation (7) and solved for a, the result is 


K. 


,1/2 


a = 


\4p t + K p y 

The equilibrium constant is related to the free energy as follows: 


( 10 ) 


l°g e K p = 


AF 

RT 


(ID 


where 


AF ° " ^ F products " ^ F reactants 


and the superscript zero designates a standard pressure of 1 atmosphere. By using the 
free energies that were calculated in the computer program for the species 


log e K p = 4.6893 
K p = 108.7772 

and with p^. = 7.1 atmospheres, from equation (10), a - 0.8905. Thus, by referring to 
equations (9), the mole numbers are 0.1095 CH 4 , 1.7810 H 2 , and 0.8905 C(s) and compare 
exactly with those determined from the computer program. This agreement indicates 
that the program is accurate. 

The second system considered was the dissociation of ammonia into nitrogen and 
hydrogen: 

NH 3 -iN 2 + |H 2 
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for which experimentally determined equilibrium compositions are given in reference 8. 
Calculated and experimental results at T = 623° K and p^. = 10 atmospheres are as 
follows : 


\ 


Mole fraction determined by — 

Species computer Equilibrium 
program constant 

NH 3 0.0772 0.0772 

N 2 .2273 .2273 

H 2 .6955 .6955 

Experiment 

0.0735 

.2282 

.6983 


The comparison, which was made at the temperature and pressure for which the greatest 
difference existed between calculated and measured results, shows that the calculations 
do not predict the equilibrium composition within the maximum experimental error of 
3 percent. (See ref. 8.) This difference apparently results from the free-energy values 
used in the calculations. Thus, as expected, the accuracy of the computer program is 
governed primarily by the free-energy input data. 

EQUILIBRIUM COMPOSITIONS FOR PHENOLIC-NYLON 

The equilibrium program was used to obtain the equilibrium compositions of the 
pyrolysis products of the low-density phenolic- nylon ablator described in reference 9. 

By the method of reference 10, the elemental composition of this material is by weight 
69.9 percent carbon, 16.3 percent oxygen, 8.12 percent hydrogen, and 5.68 percent nitro- 
gen. For the purpose of equilibrium calculations, this elemental composition completely 
specifies the material. 


Pyrolysis Productions 

Equilibrium compositions for the low-density phenolic- nylon ablator have been cal- 
culated for temperatures from 500° K to 3500° K when all the species shown in appendix C 
are considered. The results of these calculations, at each of three pressures (0.1, 1.0, 
and 10 atmospheres) are shown in figures 1(a), 1(b), and 1(c) when condensed carbon is 
included as a possible product of the thermal cracking of the pyrolysis gases. In figure 2, 
results are shown for a pressure of 1 atmosphere when the condensed species are not 
included. 

A comparison of figures 1(b) and 2 shows that the inclusion of condensed carbon has 
a very significant effect on composition. The hydrocarbons dominate the composition at 
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lower temperatures in the gaseous system of figure 2. However, when the solid carbon 
species is included, as in figure 1, these compounds do not occur until high temperatures 
are reached. 

Above 3000° K, the compositions of figure 1 become complex because of the 
increasing number of carbon hydrogen and carbon nitrogen compounds. The presence 
of carbon hydrogen compounds at these higher temperatures indicates that hydrogen is 
combining with and removing condensed carbon. Thus, carbon is removed at a much 
greater rate than would be predicted by sublimation theory alone. Therefore, the use of 
sublimation data for pure carbon to predict the mass loss rates of charring ablators at 
higher heating rates and temperatures will result in unconservative estimates. 

In the past, some investigators have attempted to compensate mathematically for 
the presence of condensed carbon by withholding some of the carbon from the original 
mixture, designating this carbon as solid or char, and proceeding with only the gas-phase 
reactions. Such an imposed constraint is somewhat arbitrary and can significantly affect 
the results. The present method eliminates the need for guessing by computing the amount 
of condensed carbon present at each equilibrium condition. 

Char density from eq uilibrium calculations.- The char density is of fundamental 
importance in ablation calculations since the surface recession rate in an oxidizing atmo- 
sphere, with no mechanical removal, is explicitly inversely proportional to the char den- 
sity at the surface. However, for quasi-steady ablation, the surface recession rate 
depends only on surface environmental conditions and free carbon density in the uncharred 
ablator, when free carbon is any carbon atoms for which no oxygen atoms exist in the 
ablator. The result assumes equilibrium between carbon and oxygen at the surface, but 
is completely independent of such factors as internal chemical reactions and char shrink- 
age. For this quasi-steady condition, equilibrium predicts the char density which will 
give the correct surface recession, at temperatures below the sublimation regime, with- 
out directly considering gas-phase reactions. Therefore, for conditions not far removed 
from quasi-steady conditions, the equilibrium approach is expected to provide a useful 
approximation, and this fact is substantiated by the experimental results of references 11 
and 12. 

By making use of this approximation, the equilibrium calculations such as shown in 
figure 1 can be used to predict char densities from the known virgin ablator densities. 

The equilibrium calculation procedure was used in conjunction with an ablation calculation 
program (ref. 13) to determine analytically the temperature distribution and char density 
of a tested ablator specimen. The specimen, which was made of low-density phenolic 
nylon (550 kg/m3), was exposed for 100 seconds to an arc-heated airstream and experi- 
enced quasi-steady ablation. The environment of this specimen was as follows: cold 
wall, nonablating heating rate, 2.26 MW/m2, enthalpy, 28 MJ/kg, and pressure, 0.03 
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atmosphere. The calculated equilibrium char density was about 318 kg/m3 for all depths 
through the char layer, since the char density is relatively constant at 59 percent of virgin 
plastic density in the temperature range of this test (2200° K surface temperature). 

The measured average density, for the entire thickness (6.3 mm) of char, was found 
to be 270 kg/m3 by air- displacement pycnometer measurements. Pyrolysis measure- 
ments of this material with thermogravimetric techniques similar to those in refer- 
ence 10 showed the residual char density to be 189 kg/m3. The residual char density 
value should correspond to the frozen-flow case where there is no cracking of the pyrol- 
ysis gases and hence no deposition of carbon by the gases to increase density. 

The pyrolysis gases are essentially frozen until they are heated to a certain tem- 
perature range. In this region, most of the gases react rapidly to equilibrium. If it is 
assumed that this region is very thin and essentially a plane, the location in the char of 
the instantaneous transition from the frozen state to equilibrium can be calculated so that 
the calculated average density would equal the measured average density. With an equi- 
librium char density of 318 kg/m3 and a frozen flow char density of 189 kg/m3, the instan- 
taneous transition was determined to be at a location about 0.7 of the char-layer thickness 
from the exterior surface. The calculated temperature at this point was 1200° K; thus, 
the transition from frozen to equilibrium conditions likely occurs at a relatively low 
temperature. 

The rate of surface recession is dependent on the char density at the surface. How- 
ever, surface recessions are usually calculated by using an average char density rather 
than the surface density, and this procedure may lead to significant errors. For example, 
the surface recession calculated for the test case just described would be 18 percent 
greater if the measured average char density were used rather than the calculated equi- 
librium char surface density. Furthermore, the overprediction of surface recession by 
using average char density could be much greater than 18 percent if the thickness of char 
layer in which chemical equilibrium occurs is small compared with the thickness over 
which nonequilibrium conditions exist. 

Pyrolysis gases.- The inclusion of solid carbon as a possible product in the pyroly- 
sis gases equilibrium has a marked effect on the mean molecular weight of these gases. 
Figure 3 shows the mean molecular weights for pyrolysis gases in chemical equilibrium. 
When solid carbon is excluded as a possible product, the mean molecular weights of the 
gases are more than twice those in which the solid phase is considered. As expected, at 
higher temperatures where the solid phase cannot exist, the two mean molecular weights 
approach each other at a given pressure. 

The aerodynamic blocking efficiency of a gas injected into the boundary layer is 
dependent on the mean molecular weight of the gas. For laminar flow, blocking is usually 
assumed to vary inversely as the 1/4 power of the molecular weight of the injected gas. 
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Thus, in the temperature range of interest (1000° K to 3000° K), the low molecular weight 
species, resulting from including solid carbon, would provide a 20-percent increase in 
blocking over that afforded by the high molecular weight gases that result when solid car- 
bon is not included as a possible product. 

The enthalpy of the pyrolysis gases in chemical equilibrium is shown in figure 4 as 
a function of temperature. The enthalpy that is plotted is the energy involved in chemical 
reactions plus the energy required to raise the gas temperatures. Thus, when moving 
from one temperature level to another, the chemical energy involved in the change in 
composition is taken into account. It is seen that at a pressure of 1 atmosphere the gases 
absorb approximately 50 percent more energy because of the chemical reactions that 
occur when solid or condensed carbon is included as a possible product. In addition, this 
energy is absorbed at a lower temperature level than when solid carbon is excluded and, 
consequently, could contribute an important increase in efficiency where surface tempera- 
tures were below 1500° K. 


The slope of the enthalpy curve is then an effective specific heat for the transpiring 
gases : 


c = 22 
P 8T 


( 12 ) 


The effective specific heat for the pyrolysis products of phenolic nylon is shown in fig- 
ure 5. There are three peaks on each curve. The one at the lower temperature is a 
result of the decomposition of methane whereas the peaks above 3000° K are a result of 
the formation of carbon hydrogen and carbon nitrogen compounds. The values of Cp in 
figure 5 can be used in an approximate numerical analysis to account for the energy 
involved in the changes in composition as well as the energy involved when the pyrolysis 
gases move from one temperature level to the next. 

There is considerable chemical energy which could have a marked effect on the tem- 
perature distribution through the char layer and on overall material performance. Calcu- 
lations to determine material performance by using a transient analysis similar to the one 
described in reference 13 have shown that heat- shield weight requirements can be affected 
by as much as 30 percent by this chemical energy (ref. 14). 


Simplified Equilibrium Determination 

An inspection of figure 1 shows that in the temperature range of 1000° K to 2700° K, 
there are four main species (nitrogen, carbon monoxide, solid carbon, and hydrogen) and 
that the proportions of these species are essentially constant over that temperature range. 
These results and other calculations for different proportions of the four elements carbon, 
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oxygen, nitrogen, and hydrogen indicated that within the temperature range of interest, the 
equilibrium composition is simply 

(1) N 2 - all the nitrogen in the four- element system 

(2) CO — in sufficient quantities to account for all oxygen in the system 

(3) CgQ^d - all carbon in the virgin plastic not appearing in carbon monoxide 

(4) H 2 - all the hydrogen in the four-element system 

With one exception, this simple composition can be assumed to apply to all materials 
which are composed of the four- element system. The exception occurs when the number 
of moles of oxygen in the material exceeds the total number of moles of carbon in the 
material. However, this situation is unlikely to exist in a practical ablation material 
because of the undesirable effect of oxygen on ablative performance. 

CONCLUDING REMARKS 

A computer program based on minimization of the Gibbs free energy has been devel- 
oped for computing multiphase equilibrium compositions of arbitrary mixtures of chemi- 
cal compounds. The program has been shown to be accurate by the comparison of com- 
puter solutions to exact solutions for simple reacting systems. Furthermore, a computer 
programing logic has been used which results in reasonable computing speed. 

The chemical equilibrium compositions of low-density phenolic-nylon ablative mate- 
rial have been computed for a wide range of temperature and pressure with the inclusion 
of solid carbon as a possible reaction product. Although not usually considered in deter- 
mining ablation pyrolysis products, the formation of solid carbon was found to have a sig- 
nificant effect on many aspects of ablation analysis. The principal effects are: 

(1) The formation of carbon- hydrogen and carbon- nitrogen compounds is delayed to 
higher temperatures. 

(2) The average molecular weight of the pyrolysis gases was much less when solid 
carbon is a reaction product. The result is that the blocking of convective heating and 
diffusion of oxygen to the heated surface is increased by as much as 20 percent for sur- 
face temperatures up to 3000° K. 

(3) The energy absorption by the pyrolysis gases due to chemical reactions is 
greater when solid carbon is included. Up to 1500° K, the energy absorption is increased 
by about 50 percent. 

(4) Subject to certain conditions, the char surface density, which is important in 
recession calculations, can be computed on a rational basis with the developed program. 
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The results also indicate that a simple equilibrium composition model (including 
four species) can be defined for the phenolic -nylon ablator and can be used in ablation 
analysis without the need of further equilibrium calculations. It is expected that similar 
equilibrium composition models could be defined for other ablative materials. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., June 18, 1969. 
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APPENDIX A 


CHEMICAL EQUILIBRIUM COMPUTER PROGRAM 
The program listing is as follows: 


PROGRAM Crl£Q ( 1 NPUT , OUTPUT , TAP 05 = F NPUT , TA PE6=QUT PUT ) 

CHEMICAL EQUILIBRIUM OF MULTIPHASE SYSTFMS BASE 0 ON THF PRINCIPLE 
OF MINIMIZATION OF THE FREE ENERGY OF THF, MIXTURE 

EXTERNAL FuFX , AFOE X 

01 MENS ION YG( 90) , YCl 10) , X! 100) , Y< 100) , SYM! 100 ) ,DELT! 100) , ENT l 100 ) , 

1 CENT! 100) ,UFrtT{ luO) ,FO< 100 ) ,3B< 10) ,CP?< 100) , F( 100) ,C! 100) , 

2 ALPHA! 10U) » XL AM ( 100 ) , F SUM ( 100) » YSUM! 100) , PMW! 100) ,WTPERC ( 100) , 

3 0SUM( 21,1 ) , IPIVUT ( 21 ) , INDEX! 2 l ,2) , COFF < 21,21 ),P(l0,10),PI(2l) 
COMMON/ BLK.1/ I , A I ( 1 00 ) , B I ( 1 00 ) , C I ( 1 00 ) , 0 I! 1 00 ) , E I l 100 ) 

COMMON/ BLR2/CP1 ( 100 ) , HI NT ( 100 ) , HO ( 1 00 ) , FHTO ( l 00 ) , FOR T < 100 ) , 

1 XMWUOO) ,SYMb< 100) ,PERC ( 100) , AA< 100, 10) , MM , ILTClvu/ 

READ INPUT uAT A - NOTE. RCAO ALL INPUTS FOR AIL GASEOUS SPECIES 

FOLLOWED bY INPUTS FOR CONDENSED SPECIES. PROGRAM WILL ACCOMMODATE 
90 GASEOUS SPLLIES AND 10 CONDENSED SPECIFS. THF NUMBER OF 
ELEMENTS CANNOT EXLhLD 10. 

5 READ! 5,2) XT, TF 1 N , T I NC , T Z ERO , CR I T 
I F ! E 0 F , 5 ) 1020,1030 

10 30 READ! 5 ,6) NG»NC , MM, NFKEQ , I DEBUG 
NN=NC, + NG 
DO 31 1=1, NN 

RT AD( 5, 8) HINT l I ) ,H0( I ) , FHTO l I ) , XMW ( T ) , Y ! I ) , SYM? { I ) 
READ(5»2)A1!I),31(I)»CI(I),0I(I),EI(I) 

31 RE AD ! 5 , 2 ) ( AA l I , J ) , J= 1 , MM ) 

READ! 9,2 )KR,P 
2 FORMAT (6E12. 4) 

6 FORMAT! 10 1 1> ) 

8 FORMAT !5F12. A, A6) 

MAXNT=50 
T = XT 
MC=NC 
MG=NG 
KX=0 

XBETA=CRI T 
500 NTEST=NFREU 
N T= 1 

DO 33 1=1, MG 
SYM! I ) =SYMb( I ) 

33 YG! I ) = Y! I ) 
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MN=MG+i 

MMN=MN 

NNN=NN 

I F { MC • EQ. 0 ) GO TO 37 
DO 34 I=MN,NN 
S YM ( I ) = SYMB ( I ) 

34 YC< I)=Y( I) 

FREE ENERGY FUNCTION 
37 00 41 I = 1 » NN 

CALL GAU$S(TZtRO,T,5G,SUM,FOFX) 

CPI ( I ) =SUM*RR 

CALL GAUSS! TZtRO.T, 50, ANS,AFOFX) 

CP2 ( I ) =ANS*KR 

DFHT ( I) = ((l./T)-C 1 ./T ZERO) ) *HI NT ( I I +(1 . /T ) *CP 1 ( I )-CP?( I ) 
FO ( I ) =FHTu ( I J+UFHT ( I ) +HO ( I )/T 

41 FORT ( I ) =Fb ( I)/(RR) 

ENTHALPY OF INITIAL MIXTURE 

ENTMI X=0. 

C M I X = 0 . 

SUMYG=0. 

DO 3 5 1 = 1, MG 

35 SUMYG=SUMYG+YG( l ) 

00 36 1=1, MG 

FNT(I) =H INT ( l ) +CP1 ( I ) 

CENT! I )=ENT( I )+HQ( I ) 

FNTMIX=ENTMIX+ENT( I ) *YS ( I) 

36 CMIX=CMIX+CENT(I)*YGm 
£ NTM I X = ENTM I X/SUMYG 

C N I X=C M I X / SUMYG 

WRITE OUTPUT FUK INITIAL MIXTURE 
WRITE! 6,204)P, T 

204 FORMAT ( I H L , 9HPRE SSURt =F 1 3. 5 , 5X , 5HT EMP= E 1 3 . 5 ) 

WRITE (6,4) 

4 F0RMAT(///&X,1HI , 1 9X, 5HF0/RT , 1 3X , 1 2HIN IT I AL Y(I)/ ) 

00 42 1=1, NN 

42 WRITE (6,14) I ,SYMii(I) .FORTCI), YU) 

14 FORMAT! IX, I o , 3X , A6 , 2fc2 0 . 8.) 
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WRITF (6fl8)ENTMlX 
WRITE (6,21)CMiX 

18 FORMAT! ///34H SENSIBLE F NTH ALP Y OF GAS M I XTUP F= E 16 . 8 ) 

21 FORMAT(/lX,3bHCHE M * SENS ENTHALPY OF GAS MIX TURE = E 1 6*81 

MASS BALANCE CONSTRAINT 

1000 DO 320 J=1,MM 
3 3 ( J ) = 0 . 0 
on 320 1=1, MG 

320 B8(J)=BB(J)+Aa( I , J ) * Y ( I ) 

IF(MC.EQ.O)GU TO 300 

DO 321 J = 1 t MM 
L = NG 

DO 321 I =MN y NN 
L = L + 1 

321 BB(J)=B8(J)+AA(L,J)*Y{ I ) 

300 Y 8 AR = 0 . 

DO 50 1=1, MG 

50 YBAR=YBAR+YG( I ) 

SET UP AND SOLVE MATRIX 

DO 10 1=1, MG 
C ( I ) =FORT ( 1 )+ALUG< P> 

F AC=YG( I )/ YbAK 
I F (FAC.LT. 1 * E-38 ) F AO= 1 • P-3 8 

10 F(I)=YG(I)*(C(I) + ALUG (FAC) ) 

IF(MC.EQ.U)GG TO 13 

J = NG 

DO LI I = Mi\ , NN 
J = J + 1 

C( I J = F 0 R T ( J ) 

11 F ( T ) = Y ( I ) *C( I ) 

1 * FFSUM = 0. 

DO A3 1=1 ,NN 
43 Ft SUM=FESUM+F ( I ) 

I F ( NT . NF . 1 ) GO TO 51 
WRITF ( 6 , 6 1 ) F t SUM 

61 FORMAT ( 1 H 27HFKEE ENERGY SUM OF M l X TUR F= F 1 6 • 8 ) 

51 DO 30 J = 1 , MM 
DO 30 K = i , MM 
ASUM=0.0 
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Dn 20 1=1, MG 

SUM= AA ( I ,J)*AA( I ,K)*YG( I ) 

20 asum=asum+sum 

30 R(J,K)=ASUM 
00 330 J=1 , MM 
AL PHA ( J ) =0 . 

DO 330 1=1, MG 

330 ALPHAt J)=ALPHA( J)+AA( I,J)*YG< I ) 
KM=MM+MC 

mm=^m+mc+i 

mmm=mm+i 

00 360 J=MMM,NM 
360 AL PH A ( J ) = 0 . 

DO 401 1=1, M M 

DO 401 J=l,l 
401 COFF ( I ,J)= ALPHA! I) 

MMC=MC +1 

IF ( MC . FQ.O ) GO TO 406 
K = NG 

DO 40? J=2 , MMC 
K = K+ 1 

DO 402 1=1, MM 
40? C OFF ( I ,J)=AA(K, I ) 

406 MCC=MMC+1 

DO 403 I =1 , MM 
DO 403 J=MCC,NM 
K=J-mmq 

403 COEF ( I * J ) =R l 1 , K ) 

DO 407 I=MMM,NM 
DO 407 J=2 , MMC 

407 COEF ( I , J ) =0. 

I F { MC . FQ . 0 ) GO TO 406 
L = NG 

00 404 I=MMH,mM 
K =0 

L = L + 1 

on 404 J=MCC,NM 
K=K+ 1 

404 COEF ( I , J ) = AA ( L , K ) 

403 DO 405 I = NM , NM 

K = 0 

DO 405 J = ,MCC,NM 
< = K + 1 


18 


■ ii ■ 


llll II 


II II II I 



o n o o o o 


APPENDIX A 


405 CDFFC I , J)=ALPrtA(K) 

DO 361 K = 1,1 
DO 361 J=i,MM 
8 SUM ( J , K ) =0 • 

DO 362 1=1, MG 

362 3 SUM ( J,K)=BSUM(J,K)+AA(I,J)*F{I) 

361 BSUM( J,K)=8SUMl J,K)+68( J) 

I F ( MC * EQ • 0 ) GO TU 409 
DO 363 K=l,l 
L = MG 

DO 363 J = MMM, KM 
L=L+l 

363 8 SUM ( J » K ) = C ( L ) 

409 DO 364 K=l,i 

DO 364 J=NM,NM 
B SUM ( J, K)=U. 

DO 364 1=1, MG 

364 B SUM ( J , K ) = B SUM ( J , K ) +F ( I ) 

CALL S IMEQC COEF ,NM, BSiJM, 1 , DET F R M , 1 P I VO T , 21 , ISCALF ) 
X Y6AR = 3SUM( 1,1) 

I F ( MC • EQ • 0 ) GO TU 103 

NEW X-S FOR CONDENSED SPFCIFS 


J = MG 

DO 101 1=2, MMC 
J=J + 1 

101 X ( J ) = B SUM { 1,1) 

NEW X-S FOR GASEOUS SPtCIES 
103 J = 0 

DO 102 I = MCC , NM 

J = J+l 

102 PI ( J ) = BSU M ( 1,1) 

00 60 1=1, MG 

60 FSUM1I )=-F(I )+YG( I >*XYBAR 
DO 110 1=1, MG 
P SIJM= 0 • 

00 120 J= 1 , MM 

120 P SUM = PSUM + P I (J)*AA(I,J) 

Y S UM { I )=PoUM*YG( I ) 

110 X ( I )=FSUM( I J+YSUMU) 
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C 

C LAMBDA AND DIfttCT 10NAL DERIVATIVE 
C 

XLAMBD = 1. 

310 DO 86 1=1, NN 

DFLT ( I ) =X { I >-Y( I ) 

I F ( IDEBUG.NE.OJPRINT 1005,DELT(I) 

1005 FORMAT! lX,tl5.5) 

IFIOELT(I) .Gt.O. )GO TO 36 
XL AM ( I ) = Y ( I ) / ( -DELT ( I ) ) 

XLAMRD=AMIN1 { XLAMBD, XL AMI I ) ) 

86 CONTINUE 
DEBAR=0. 

00 87 1=1, MG 

87 DEB AR = DE BAR+DLL1 I I ) 

93 DFDL = 0 . 

DO 88 1=1, MG 

96 FAC=(YG( I ) +XLAMBD*DELT ( I ) ) / I Y8AR+XL AM BD* OF BAR ) 

IF! FAC.GT.O. )G0 TO 82 
XLAM3D=.5>S'XLAM6D 
GO TO 96 

82 DFDL = DFOL-t-DLLT (I) * ( C ( I ) + ALOGI ( YG ( I ) +XL AMSD*DE L T ( I ) ) / 
1 (YRAR+XLAMGU*OeGAR) I ) 

38 CONTINUE 

1 F ( MC . E 3. 0 ) GU TO 84 
00 83 I = Mi«i,NN 

83 OFDL = OFDL + C( l )*OcLT ( I ) 

84 IFIOFOL.Lr. 0.)Gu TO 69 

I F ( XLAMBD.GT.O.E— 38)GO TO 91 

KX=KX+1 

GO TO 94 

91 XLAMBD=.5*XLaMBD 
GO TO 93 

89 IF! IDEBUG.Nt.OlPRINT 1021,XIAM80 
L021 FORMAT! //IX, /HLAMBDA=E l 5 .5// ) 

C 

C NEW Y-S FOR GASfcUUS AND CONOCNSFD SPECIFS 
C 

94 00 76 1=1, NN 

Y ( I ) =Y< I H-XLAM60*DfcLTII) 

I F ( IDEBUG.Nt.OlPRINT 1031,SYM! I ) , Y ( I) 

1031 FORMAT ( IX, A6, LI. 6.5) 

76 CONTINUF 
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DO 77 I = 1 f MG 
77 YG ( f >=Y< I) 

I F { MC • EQ. 0 ) GO TO 76 
J -NG 

DO 74 I=MN,NN 
J = J + l 

74 YC( J)=Y{ I) 

MOLE PERCENT 

7 B SUMY-O. 

DO 370 1=1, NN 
370 SUMY-SUMY+ Y ( I ) 
on 340 1=1, NN 

340 PERC U»=( Y( I)/SUMY)*100. 

1 F ( MC . EQ . 0 ) GO TO 351 
J = NG 

DO 350 I = MN , NN 
J=J + 1 

350 P EKC ( J ) =PLKC ( I ) 

TEST FOR CONVERGENCE 

351 8ETA= 0 . 

DO S 5 1 = 1, NN 

85 BETA = B6TA + ABStUELT( 1)) 

IFCBETA.LT .XBETAJGU TO 800 
IF(KX.GT.O)GL) TO 560 
l F ( NT .FQ.NTEST I GO TU 570 
GO TO 560 

TEST FOR SIGNiFICANLE OF SPECIES 

570 NTFST-NTFST+NFR EQ 
590 IHMG.FQ.OJGU TU 48 

CALL TEST( YG, 1 , MG , NMG ) 

48 l F ( MC • EO • 0 ) GO TO 47 

CALL TEST ( YL , M MN , NNN , MNN ) 
MC=MNN-NG 
47 MG=NMG 
MN=MG+1 
NN-MC+MG 
NNN=NG+MC 
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00 561 1=1, MG 
S YM ( I ) =SYMB( 1 ) 

561 Ym=YG(I) 

1 F t MC • EQ. 0 ) GU TO 5t>0 
J=NG 

00 562 I = MN» NN 

J=J+l 

S YM ( I ) =SY MB ( J ) 

562 Y( I ) = YC ( J ) 

560 I F ( NT.GE .MAXNT )6U TO 600 
NT =NT+ 1 

IFJKX.GT.OJGU TO 300 
GO TO 1000 

600 XBFTA=XriETA+CRlT 
MAXNT= MAXNT +50 
N T = NT + 1 
GU TO 100U 

800 NN=MC+NG 
«X=KX+1 
MG=NG 
MN=MG+ 1 

00 563 1=1, MG 
S YM ( I )=SYMd( I ) 

563 Y( I ) =YG ( I ) 

I F ( MC . FO.OJUU TO 56 J 
00 564 I = MN , NN 
S YM( I ) =SYM3( I ) 

564 Y(I ) =YC ( I ) 

567 I P ( KX .LT.2 )NT =NT+1 

IF(KX.LT.2)Gu TO 1000 

WEIGHT PERCENT 

801 T PMW= 0 • 

DO 30 1=1, NN 

PMW ( I ) =P p KC ( I ) *Xrtw ( I ) 

80 TRMW=T PMW+PMw ( I ) 

00 70 1=1, NN 

7 0 WTPERC U >=tPMw( I ) / TP MW ) * 100. 
KKIB=NN+1 

ENTHALPY OF ti*U I L I BR IUM MIXTURE 
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CMI X=0. 

ENTMI X=0. 

XSUMY=0. 

00 580 1=1, MG 
580 XSUMY = XSUMY+-YG( I ) 

00 73 1=1 ,NN 
F NT ( I ) =H I NT ( LIHPltl) 

73 CFNTl I )=ENT( I >+H0( l ) 

DU 75 1=1, MG 

ENTMIX=ENTM1X+ENT( 1 )*YG< I ) 

75 CMIX=CMIX+CENI ( 1 )*YG( I ) 

ENTMI X=FNTMI X/XSUMY 
C M I X=C M I X/XSUMY 

MOLECULAR WEIGHT OF EQUILIBRIUM MIXTURF 

XMWMI X=0. 

00 81 1=1, MG 

81 XMWMlX=XMwMlX+XMW( 1 )*YG( I ) 

XMWMI X=XMwMl X/XSUMY 

WRITE OUTPUT FOR EQUILIBRIUM MIXTURF 

WRITF(6,lo)NT 

16 FORMAT! / 1 X , 2 7HNT= NO. ITERATIONS REQUIRFD=I5/ 1 
WRI TE ( 6,2ul) XbETA 

201 FORMAT! / 1 X , 5HbETA=E12 .4/ ) 

WRITE (6,29) 

29 FORMAT ( //47X , 4HM0LE ,17X,6H WEIGHT) 

WR I TE ( 6 , 32 ) 

3? FORMAT ( 6X , 1H I ,21X,4HY( I } , 14X, 7 H PERCENT , 14X,7HPERCFNT / ) 

DO 44 1=1 , NN 

44 WRITE (6,19 ) 1 , SYMb ( I ) , Y ( I ) , PFRC ( I ) , WTPFRC ( T > 

19 FORMAT! IX, lt>,3X»A6,3E?O.R) 

WRITE (6,12) 

12 FORMAT! /// 12X, 17HSENSIBLF ENTHALPY , 3X , 20HCHFM + SENS ENTHALPY/) 
DO 79 1=1 ,NN 

79 WRITE (6,24 1 5YMBII) , ENT ( I ) .CENT ( I ) 

24 FORMAT ( IX , A6.2L2G.8) 

WRITE (6,18 ) ENTM1X 
WRITE (6,21 ) CM 1 X 
WRITE (6, 17) XMWMI X 

17 FORMAT ( /IX, 32H MOLECULAR WEIGHT OF GAS M I XTURE =E 1 6 . 8/ / ) 
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WR I TE { 6 f 6 1 )F ESUM 

I NCREMFNT TtM PERATUR t AND START NEW CASE 

IF (T.LE.TFINiGU TO 5 
T-T-TINC 
X BET A=CR I T 
K X=0 
MC=NC 
NN= NC + NG 

I F ( MC • EQ • 0 ) GO TO 500 
DO 52 I =KKL B » NN 

52 YC(I)=l.F-03 
DO 53 I=KKLb » NN 

53 YU ) = YC ( I ) 

GO TO 500 

1020 STOP 
END 


FUNCTION FUFX(TEMP) 

CGMMON/BLM/I t AI (100) t B I ( 1 00 ) , C I ( 1 00 ) , D I ( 1 00 ) f E I ( 1 00 ) 

FOFX=AI ( I)+di ( I)*TEMP+CI ( I )*TEMP**2+DI ( I )*TEMP**3+EI ( I )*TEMP**4 

RETURN 

END 


FUNCTION APUFX(TEMP) 

COMMON /BLK1/I »Al( 100) , «I < 100) ,CI < 100), DI< 100) ,EI < 100) 

A FOE X = ( AT U)+bI { 1 ) *TEMP+CI ( I )*TEMP**2 + DI ( I )*TEMP**3+F Id) *T EMP**4) 
1/TFMP 
RETURN 
END 
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SUBROUTINE: TEST ( V AR » I N» N , M ) 

COMMON/BLKI/I , A l ( 100) ,BI(100I,CI(100),DI(100),EI(100) 
C0MM0N/BLK2/CPU 100) ,HINT( 100 ) , HOI 100) ,FHTO( 100) , FORT ( 1^0) , 

1 XMW ( 100) , SYMBI 100 ) ,PERC (100 ) , AA« 100, 10) ,MM, i;i £LT(1 0) 

DIMENSION VAR (100) , Y X( 100) , X YC 1 00 ) , XA I ( 1 00 ) , X B I ( 1 00 ) , XCI ( 1 00 ) , 

1 XDI ( 109) ,XEl ( 100) tXCPIC 100) ,XHINT( 100 ) , XHO ( 1 00 ) , XFHTO ( 1 00 ) , 

2 XFORT (100) , XXMWI 100) ,XSYMB( 100 ) , X P FR C ( 1 00 ) , X AA ( 1 00 , 1 0 ) 
C0RR=l.E-50 

00 531 I=IN,N 
IF <VAR( I ) ,NE.O.)GO TO 531 
VAR(l) =COKR 
C QRR=CORR+ 1 . E— 50 

531 CONTINUE 

00 540 J= I N » N 
YMAX=0. 

DO 520 L= I N , N 

520 YMAX=AMAX1(YMAX, VAR ( L ) ) 

90 530 K= I N » N 
YX(K)=VAR(K.)/YMAX 

IF(ABS(YX(K)-1. ).GT.1.F-7)G0 TO 530 
V AR ( K ) =0. 

X Y ( J )=YMAX 
XAI ( J > =AI ( K) 

XBI ( J )=BI ( K ) 

XCI ( J ) =C I ( K ) 

XDI ( J ) =DI ( K) 

XEI ( J ) =F I ! K.) 

XCP1 ( J ) =CP 1 ( K) 

XHINTl J)=H1NT(K) 

XHO ( J ) = H0 ( K. ) 

XFHTO( J)=FHTU(K) 

XFORT ( J ) =FURT { K ) 

XXMW( J) =XMw( K) 

XSYMBt J) =SYMB( K) 

XPERCt J)=PERC(K) 

lt(j LTC,) 

00 532 1=1, MM 

532 XAA( J , I )=AA( K, I ) 

GO TO 540 

530 CONTINUF 
540 CONTINUF 

00 510 I = 1 N , N 
A I ( I ) =X A I ( l ) 
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. lv( i )=;: /:lk i ) 

b i ( k=xbi ( n 

ci ( I)=XCI ( I) 

01(1 l = XDI ( I) 

E I ( I ) =XE I ( I ) 

S YMB ( I )=XSYMb( I ) 

FORK I ) =XFUk T ( I ) 

CP1U ) =XCP 1 ( I) 

H I NT ( I ) =XH 1 NT ( I) 

HC< I ) = XHO ( I ) 

F HTO ( I ) = XFHTO( [ ) 

XMW( I ) =XXMw ( 1 I 
PERC ( I )=XPeRb ( I ) 

00 590 J=1 , MM 
590 A A ( I , J )=XAA{ I , J ) 

510 V AR ( I ) = X Y ( I i 
M = N 

DO 550 I=IN,N 

IF(PFRC( I ) .iiT. I.t -5 ( I ) . sT.o . )•■ j Tv. 550 

M=M - 1 

550 CONTINUE 
RETURN 
END 


SUBROUTINE GAUSS ( A , B , N , SUM ,F0FX ) 

DIMENSION U l 5 ) * R ( 5 I 

U(l) =.4255o2830509184 

U( 2> =.253302302985376 

U( 3) =. 160295215550488 

U ( 4)=.06746o3l6655508 

U( 5)=. 013046735741414 

R ( 1)=. 1477o21 12357376 

R(2)=. 1 34633359654998 

R ( 3) = . 109543181257991 

R (4) =.074725o 745 75290 

R (5)=. 033335672154344 

SUM=0. 0 

IF(A.PQ.B) RETURN 
F I NF = N 

OEL TA=F INL/ ( d-A) 

DO 3 K=l,N 
X I = K— 1 

F I NF = A +X I /DELTA 

00 2 11= 1,5 

UU=U( I I ) /UlLTA + FINE 

2 SUM=R( II )*HJFX< UUI + SUM 
00 3 1=1,5 

UU=( 1 . 0-U (L) I/DELTA +F INF 

3 SUM=R (L )*FUFXI UU1+SUM 
SUM= SUM/ DELTA 
RFTURN 

END 
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C SOLUTION OF SIMULTANEOUS LINEAR EQUATIONS 

********* DOCUMENT date 08—0 1 —66 SUBROUTINE REVISED 08-01-08 ******* 
SUBROUT INE S I MCQ ( A * N * b * M * DE TERM * I P I VOT * NM AX * I SCALE ) 

C 

D I MENS ION IPI VOT ( N ) * A ( NMAXtN ) * B (NMAXiM ) 

EQUIVALENCE ( IROWt JROW) ^ ( I COLUM * JCOLUM ) * ( AM AX * T * S WAP ) 

I N I T I AL I ZAT I ON 
C 

5 ISCALE=0 

6 R1 =1 0*0**100 

7 R2 = 1 .O/Rl 

10 DETERM= 1 • 0 
15 DO 20 J=1 *N 
20 IPIV0T( J)==0 
30 DO 550 1=1 tN 

C 

C SEARCH FOR PIVOT ELEMENT 

C 

40 AM AX = 0 • O 
45 DO 105 J=UN 

50 IF ( IPI VOT ( J )-l ) 60 <105*00 

60 DO 100 K= 1 <N 

70 IF (IPIVOT(K)-l ) 80<100<740 

80 IF (ABS ( AMAX )~A°S ( A ( J*K ) ) ) 85*100*100 

85 1R0W®J 

90 I COLUM = K 

9b AMAX=A(J*K) 

100 CONTINUE 
10b CONTINUE 

IF (AMAX) 1 10* 106* 1 10 
106 DETERM = 0 • O 
I SCALE^O 
GO TO 740 

110 IPI VOT ( I COLUM ) = I P I VOT ( I COLUM ) + 1 
C 

C INTERCHANGE ROWS TO PUT PIVOT ELEMENT ON DIAGONAL 

C 

130 IF ( IROW-ICOLUM) 140,260*140 
1 4 0 DE TERM = — DETERM 
1 50 DO 200 L= 1 *N 
160 SWAP = A ( I ROW * L ) 

170 A ( IROW*L >=A ( I COLUM *L ) 

200 A ( I COLUM * L ) = S W Ap 
20b IF (M> 260*260*210 

210 DO 250 L= 1 * M 
220 SWAP = B ( I ROW* L ) 

23U B ( I ROW * L )=B ( I COLUM, L) 

250 B ( ICOLUM*L )=SWAP 

260 P I VOT = A ( I COLUM * I COLUM ) 

IF (PIVOT) 1000*106*1000 
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scale the determinant 

10C0 P I VOT I =PI VOT 

10 05 IF (ABS(OETtiRM)-Rl )10 30 *1010*1010 
1010 DETERM = DETERM/R 1 
ISCALE= ISCALE+1 

IF ( AB S ( DET ~ RM ) -R 1 ) 1060* 1020* 1020 
1020 DETERM = DETERM/R 1 
!SCALE=ISCALE+1 
GO TO 1060 

10 30 IF (ABS CD'S TERM)- R2 ) 104 0* 104 0. 10 60 
104 0 D E T t R M = DE Tc. RM* R 1 
1 SCALE= iscale-i 

IF ( ABS ( UET-RM ) -R2 ) 1 050 * 1050*1060 
1 0 50 DtTERM = Dc.TLRM*Rl 

iscale= iscalf:-i 

1 060 I F ( ABS ( P I VOT I ) -R 1 ) 1 090 * 1 0 70 * 1 070 
107' P I VOT I = P I VOT I/Rl 
I SCALE = lSCALt+1 

I F ( ABS (P 1 VOT I ) -R1 >320,105 0*108 0 
1080 P I VOT I =P I VOT I /Rl 
1 SCALE = I SCALE + 1 
GO TO 320 

1090 IF ( ABS ( P I VOT I ) -R^ ) 2000 * 2000 * 320 
2000 PlVOTl=PlVOTI*Rl 
ISCALE= ISCALE-1 

I F ( ABS (PIVOT 1 )-R2 >2010*20 10*32 0 
20 I 0 P I VOT 1 = P I VOT I *R 1 
ISCALE= ISCALE-l 
32C DE7ERM = DETtiRM*P I VOT I 
C 

C DIVIDE PIVOT ROW BY PIVOT ELEMENT 

C 

34 G DO 351 L = 1 * N 

3^1 IF ( I P I VOT (L >- 1 ) 350 ,351 *740 

350 A { I COLUM * L ) = A ( I COLUM * L > /P I VOT 

351 CONTINUE 

355 IF (M) 380*380*360 

360 DO 370 L = 1 * M 

370 B { I COLUM, L )=t>( I COLUM *L ) /PI VOT 
C 

C REDUCE NON-PIVOT ROWS 

C 

380 CO 550 LI = 1 * N 

390 IF (L 1 - I COLUM ) 400*^50*400 

4 C 0 T = A (L 1 * I COLUM ) 

430 DO 451 L = 1 * N 

431 IF ( IPI VOT (L >-l > 450*451*740 

450 A f L 3 *L ) =A ( LI ,L ) -A ( I COLUM *L ) *T 

451 CONTINUE 

455 IF (M) 550*550*460 

460 DO 500 L = 1 * M 

500 B ( L 1 * L ) =b (L 1 *L)-b(I COLUM * L ) * T 
55 p CONTINUE 
740 RETURN 
END 
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The 

XT 

TFIN 

TINC 

TZERO 

CRIT 

NG 

NC 

MM 

NFREQ 

IDEBUG 

HINT (I) 

HO(I) 

FHTO(I) 

XMW(I) 

Y(I) 

SYMB(I) 


INPUTS FOR EQUILIBRIUM PROGRAM 
inputs for an equilibrium program are as follows : 
initial temperature at which equilibrium is to be calculated, °K 
final temperature at which equilibrium is to be calculated, °K (XT = TFIN) 
temperature increment at which equilibrium is to be calculated, °K 
298.2° K 

N 

convergence criterion ^ |xj - y||< CRIT 

i=l 

number of gaseous species 

number of condensed species 

number of elements 

frequency of test for negligible species 

trigger for printing debug information; IDEBUG not equal to 0 prints this 
information 

sensible enthalpy for standard state of ith species, cal/mole 

chemical energy at 0° K for standard state of ith species, cal/mole 

sensible free energy for standard state of ith species at 298.2° K, cal/mole 

molecular weight of ith species 

initial guess for the mole number of ith species 

alphameric name of ith species 
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AI(I), BI(I), CI(I),''* coefficients of equation for specific heat of ith species. The equa- 
DI(I), EI(I) J tion has the forna: c p = A + BT + CT2 + DT3 + ET4 

AA(I,J) formula number indicating number of atoms of elements j in a molecule of 
species i 

RR universal gas constant, 1.987 cal/mole-°K 

P pressure, atmospheres 
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LIST OF SPECIES 

A list of species for which data cards have been prepared is presented herein. 
Data for the following species are contained in reference 15: 


Gases 

Gases 

Solids 

C 

N 

C 

c 2 

n 2 

Si solid 

C3 

NH 

Si liquid 

CH 

nh 3 


ch 2 

NO 


ch 3 

no 2 


ch 4 

n 2 o 2 


c 2 h 2 

0 


c 2 h 4 

o 2 


CN 

OH 


c 2 N 2 

Si 


CO 

Si 2 


co 2 

Si 3 


H 

SiH 


h 2 

sm 4 


HCN 

SiN 


HCO 

SiO 


H 2 o 

Si0 2 



The species from reference 16 are: 


c 2 h 6 

c 3 h 8 

CHO 

CH 3 OH 

c 3°2 

COH 2 

h 2 o 2 

°3 
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The following species were obtained from reference 17: 

C 2 H 

c 3 h 

c 4 h 

c 4 h 2 

The species from reference 18 (solids) are: 

SiC 

Si0 2 
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Figure 2.- Equilibrium composition when only the gaseous phase is allowed. Pressure, 1 atmosphere. 
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